function [ S ] = CubicSplineInterpolation(fx, points)
f = inline(fx);
Y = f(points);
N = size(points, 2) - 1;
init = zeros(1, N);

S = struct('xp', num2cell(points(1, 1 : N)), 'xi', num2cell(points(1, 2 : N + 1)), ...
    'a', num2cell(init), 'b', num2cell(init), 'c', num2cell(init), ...
    'd', num2cell(init));

A = zeros(N);
B = zeros(N, 1);
A(1, 1) = 2/3 * (S(1).xi - S(1).xp) + 2/3 * (S(2).xi - S(2).xp);
A(1, 2) = (S(2).xi - S(2).xp) / 3;
B(1) = (Y(3) - Y(2)) / ((S(2).xi - S(2).xp)) + (Y(2) - Y(1)) / (S(1).xi - S(1).xp);
for i = 2 : N - 1
    hi = S(i).xi - S(i).xp;
    hi_1 = S(i + 1).xi - S(i + 1).xp;
    
    A(i, i - 1) = hi / 3;
    A(i, i) = 2/3 * hi + 2/3 * hi_1;
    A(i, i + 1) = hi_1 / 3;
    B(i) = (Y(i + 2) - Y(i + 1)) / hi_1 + (Y(i + 1) - Y(i)) / hi;
end;
A(N, N - 1) = (S(N - 1).xi - S(N - 1).xp) / 3;
A(N, N) = 2/3 * (S(N - 1).xi - S(N - 1).xp) + 2/3 * (S(N).xi - S(N).xp);
B(N) = (Y(N + 1) - Y(N)) / ((S(N).xi - S(N).xp)) + (Y(N) - Y(N - 1)) / (S(N - 1).xi - S(N - 1).xp);

c = ThomasAlgorithm(A, B);

for i = 1 : N 
    S(i).a = Y(i);
    S(i).c = c(i);
    hi = S(i).xi - S(i).xp;
    if (i == N)
        ci_1 = 0;
    else
        ci_1 = c(i + 1);
    end
    S(i).b = (Y(i + 1) - Y(i)) / hi - c(i) * hi - ((ci_1 - c(i)) * hi) / 3;
    S(i).d = (ci_1 - c(i)) / (3 * hi);
end;
end

